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ABSTRACT 


The solution of small-strain elastic-plastic stress analysis problems by the p- 
version of the finite element method is discussed. The formulation is based on 
the deformation theory of plasticity and the displacement method. Practical real- 
ization of controlling discretization errors for elastic-plastic problems is the main 
focus of the paper. Numerical examples, which include comparisons between the 
deformation and incremental theories of plasticity under tight control of discretiza- 
tion errors, are presented. 
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INTRODUCTION 

This paper is concerned with application of the p-version of the finite element 
method to elastic-plastic stress analysis problems with emphasis on the deforma- 
tion theory of plasticity. Our interest in this subject is motivated by the following 
considerations: 

(1) The effects of a single overload event on structures made of ductile materials 
are of substantial practical importance. Such effects can be well represented 
by mathematical models based on the deformation theory plasticity [l], [2]. 
This is illustrated by four examples. 

(2) The propagation of cracks in strain-hardening materials is generally correlated 
with the J-integral. The J-integral is based on the deformation theory of 
plasticity [3], [4]. 

(3) The p-version is not susceptible to Poisson ratio locking and hence correct 
limit loads are obtained. In the conventional (h-version) locking occurs when 
the displacement formulation is is used. For this reason alternative formula- 
tions, generally known as mixed methods, must be employed. See, for exam- 
ple, [5]). 

(4) Realistic mathematical models of real physical systems must have a capability 
to provide initial estimates for the effects of nonlinearities at a low computa- 
tional cost. The deformation theory of plasticity serves this purpose well. 

(5) Adaptive control of discretization errors is more important in the case of 
nonlinear problems than in the case of linear problems because the initial 
discretization may not be adequate throughout the solution process, hence 
errors may accumulate in the course of iteration. The p-version, which uti- 
lizes hierarchic finite element spaces, is well suited for controlling discretiza- 
tion errors: The number of degrees of freedom can be increased substantially 
without mesh refinement. It is particularly advantageous to use adaptive 
p-distributions in nonlinear cycles because the gains in performance, as com- 
pared with unadapted schemes, are multiplied by the number of iteration 
steps. 



(6) From the point of view of implementation, the data storage requirements for 

the deformation theory are much smaller than for the incremental theory. 

The p-version of the finite element method became well established during 
the 1980’s. Its theoretical basis is now thoroughly developed and its performance 
characteristics are extensively documented. We refer to [6] and the references 
listed therein. With the exception of adaptive hp-extensions in fluid dynamics 
(see, for example, [7]), virtually all documented applications of the p-version have 
been to linear problems, and particularly to problems belonging to the follow ing 
two categories: 

Category A: The exact solution is analytic on the entire solution domain and its 
boundaries. 

Category B: The exact solution is analytic on the entire solution domain and its 
boundaries, with the exception of a finite number of points (in three dimensions 
finite number of points and lines). The points where the solution is not analytic 
are called singular points. 

The p-version is effective for problems in Categories A and B because expo- 
nential convergence rates can be achieved, within the range of accuracy normally 
expected in engineering practice, with simple finite element meshes. 

The effectiveness of the p-version depends on the smoothness of the underlying 
exact solution Hex and the design of the finite element mesh with respect to Hex- 
Certain types of nonlinearities, such as material nonlinearities associated with 
nonlinear elasticity and the deformation theory of plasticity, and even a broad 
class of problems solved by the incremental theory of plasticity, do not perturb 
the smoothness of the underlying exact solution significantly. Therefore the p- 
version is an effective method for solving such problems. In fact, the performance 
characteristics of the p-version can be expected to be substantially the same as 
in the case of linear problems belonging in categories A and B. This expected 
behavior has been confirmed with respect to a set of benchmark problems, four of 
which are presented in this paper. 

The assumptions on which the deformation theory of plasticity is based; the 
elastic-plastic compliance matrices for the cases of plane stress, plane strain and 
axisymmetric problems; an outline of the algorithmic procedure and examples are 
presented in this paper. 
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FORMULATION OF THE MATHEMATICAL PROBLEM 
The cases of plane stress, plane strain and axially symmetric problems, that is, 
two-dimensional formulations, are considered in the following. The formulations 
are based on the displacement method. 

Notation. 

The components of the displacement vector u are denoted by u x = u*(x,y) and 
u v = «y(*,y) The components of the small strain tensor, by definition, are: 



In addition, 7 * v = 2e Itf will be used to represent the usual engineering definition 
of strain. The elastic (resp. plastic) strains will be indicated by the superscript 
e (resp. p). The three principal strains are denoted by «i, e 3 , e 3 . The equivalent 
elastic strain is defined by: 

«* " $ " £ 5) 3 + ( £ S " £ i) 2 + (4 - £ !) 2 (2a) 

where v is Poisson’s ratio. The equivalent plastic strain is defined by 

* d = + M 

and the total equivalent strain is, by definition, 

? d A f ?« + **>. ( 2c ) 

The uniaxial strain at the onset of yielding is denoted by ty. 

The stress tensor components are denoted by o-*, <r y , a, , r xy . The three prin- 
cipal stress components are denoted by <n, a 3 , a 3 . The equivalent stress is defined 
by: 

9 VI 0 '! — ^a) 3 + (a 3 — <r 3 ) 2 + (<r 3 — <ri) 2 (3) 
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The components of the stress deviator tensor are denoted by a x , a y , & M , By 
definition, 


cr, =a x - -{<r, +<r y + a M ) 

1 , 

<Ty =<Ty -j{<T M +<Ty+ O ,) 

1 , 

<r, =cr, - -(<7, + tr„ + <r M ) 

T X y —T X y 


(4a) 

(46) 

(4c) 

M 


The second invariant of the stress deviator tensor is denoted by J a and is 
defined by: 

J2 = \(al + 5l + a 3 M ) + ^y. (5) 

In the case of axial symmetry the independent variables are denoted by r, 6, 
z instead of x, y, z. In the one-dimensional case (i.e., uniaxial stress state) the 
subscripts are omitted. 


Assumptions. 

The assumptions on which the formulation of the mathematical problem is 
based are described in the following. 

Assumption 1: 

The strain components are much smaller than unity on the solution domain 
and its boundary, and the deformations are small in the sense that equilibrium 
equations written for the undeformed configuration are essentially the same as the 
equilibrium equations written for the deformed configuration. 


Assumption 2: 

The total strain is the sum of the elastic strain and the plastic strain. 
Referring to Fig. 1, in the case of uniaxial stress state the stress-strain law is: 



where E, is the secant modulus. Since the elastic part of the strain is related to 
the stress by Hooke’s law: 



- 4 - 



where E is the modulus of elasticity, we have: 


€* = 



( 6 ) 



Fig. 1. Typical uniaxial stress-strain curve. 


Assumption 3: 

The absolute values of the stress tensor components are non-decreasing and 

the stress tensor components remain in a fixed proportion as the deformation 
progresses. 

Assumption 4: 

The plastic strain tensor is proportional to the stress deviator tensor. 
Assumptions 3 and 4 allow generalization of the uniaxial stress state for which 
experimental information is available to two and three dimensions. In the case of 
uniaxial stress state a = 2<r/3 and hence eq. (6) can be written as: 



In two-dimensional problems: 




r V* ’ 


3 / 1 1 \ 

By 





2 \E. EJ ' 

Oz 



^f X y - 


(7a) 


( 76 ) 
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Remark: 


In the incremental theory of plasticity based on the von Mises yield criterion 
the following assumption, analogous to Assumption 4, is made: Increments of the 
tensor components of the plastic strain are proportional to the first derivatives of 
J 3 with respect to the corresponding components of the stress tensor. The first 
derivatives of with respect cr,, <r„, <r, and r xy can be shown to be equal to 
5, and respectively. Thus, for example; 


= dX = dXa a . 

0(7 x 


In the one-dimensional case: 


dt p = d? = d\a = -d\9 
3 


hence: 


, p 


(>) 


Analogous relationships hold for the plastic increments of each component of the 
strain tensor. 


The elastic-plastic material compliance matrix in the case of plane stress. 

Using the definition of the stress deviator, given by (3), and the relationship 
between the plastic strain and deviatoric stress (7b), we have: 


f ‘S ' 


■ 2/3 

-1/3 

0* 


' ) 

< P V 

3 / 1 l\ 

' 2 \E. Ej 

-1/3 

2/3 

0 

i 

<Jy 

' Ixy * 


. 0 

0 

2. 


i T X y ) 


Using {e} = {£*}+{e p } a relationship is obtained between the total strain components 
and the stress tensor: 


[ ^ ) 


' 1 -v 0 


• 1 -1/2 O' 

\ 




i 


E-E. 





1 W 1 

E 

-v 1 0 

+ * 
E,E 

-1/2 1 0 


i 

CTy 

1 i 


. 0 0 2(1 + v ) . 


. 0 0 3. 

) 


< Txy ' 



The matrix in the brackets is the elastic-plastic material compliance matrix which 
is readily invertible to obtain the material stiffness matrix. 

Remark: 

Referring to Fig. 1 and the definitions for the equivalent plastic strain and 
equivalent stress, it can be easily shown that 


( E-EJ €»> P 
EE, ~ cr ~ & 


In view of Eq. (8), in the incremental theory of plasticity the equation analogous 
to (9a) is: 


( de x ' 


" 1 

-v 0 ' 

[ 

1 

' “ E 

— 1/ 

1 0 

l d'fxv , 


. 0 

0 2(l + i/). 


do* \ 

• 1 

-1/2 

O' 


’ <T. ' 







M + - 

-1/2 

1 

0 

i 

<Ty 

dr xv ) 

. 0 

0 

3. 


k T X y 


(96) 


The elastic-plastic material compliance matrix in the case of plane strain. 
In the case of plane strain we have 

e , = e ; + e p = o 


where 


Therefore: 


and 


e * = E “ va * ~ 




l 

° x ~ 2 


(*»+**) 


= a* - J (<T« + Vy + 0*) 


- + *v + \ [l - ^-(1 - 2v)] [<r M + <T y ) j 

\ + if i 1 - H '■-[?- if < i - H *v 
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Similarly: 


*• " ■ [\ - ST* 1 - 2v) ] + [l + ?T (1 ■ 2v) \ "■'* 

The elastic strain components in terms of the stress components are: 

= i( l ~ * + sf (i - 2 " ) ) - i(!-H , i - 2 ‘' ) ) < '’ 


ixy 


2(1 +id r 

~ E - 


xy 


and the plastic strain components in terms of the stress components are: 

_ E-E. \3 IE,. A E — E, f 3 IE, _ A 

** " E.E [4 + 4 E ^ 2v A ° * [4 4 E ^ 2 ^J 

E-E, \ 3 lE,. t „ .1 £-£, fS , l£ l(< . A 

e *~ E.E [4 4E^ 2 ^J <T » + E.E [4 4 2? * 


E-E. 

4, •*-=-=*■* 


(10a) 

(106) 

(10c) 


E.E 


xy 


( 11 °) 

(116) 

(11c) 


Combining equations (I0a,b,c) and (11a, b,c), the elastic-plastic material compli- 
ance matrix can be written in the form: 


[C] d = 


Cu 

C12 

0 


C 12 

C X1 

0 


0 

0 

C 33 


where: 



E-E. 

j i 

EE. 


Nf(-H 


Cl2 = ~l (f“H {1 ~ 2, ' ) ) 


E-E, 

EE. 



C33 = 


2(1 + 1/) t . E-E, 
E EE. 


To obtain the elastic-plastic material stiffness matrix, (C| is inverted. 



The elastic-plaatie material compliance matrix in the case of axial symmetry. 

In the axisymmetric case the elastic part of the the radial, circumferential 
and axial strain components are related to the corresponding stress components 
by Hooke’s law: 

« 1 t i 

c r = — (<r r - va 9 - v<t m ) 

« 1 , 

e, = — {-u<r r + <r 9 - v<t m ) 

€ 'm = \ (~ va r - V<T 9 + <T,) . 

The plastic strain components are related to the stress by: 


ren 

■ 2/3 

-1/3 

-1/3- 


r Vt ' 

9 | 2 \E. Ej 

-1/3 

2/3 

-1/3 

i 

<r 9 

uj 

.-1/3 

-1/3 

2/3 . 


>. j 


The elastic-plastic compliance matrix is of the form: 


[CJ 


’Cu Cm C\2 
Ci a Cn Cm 
Cm Cm Cm 
. 0 0 0 


0 ' 

0 

0 


C44 J 


where: 



Cm = — 


2 EE, 


(E — E,( 1 — 2i/)), C 4 . 


2(1 ±v) , E-_E. 

E EE. 


Since the elastic-plastic compliance matrix has a special structure, its inverse can 
be readily computed to obtain the elastic-plastic stiffness matrix. 



OUTLINE OF THE SOLUTION ALGORITHM 
Following is an outline of the procedure used in solving the elastic-plastic prob- 
lems based on the deformation theory of plasticity described in the next section. 
The procedure is known as direct iteration. The iteration number is represented 
by a superscript in brackets. 

1. Obtain a linear solution. Ensure that the relative error in energy norm is 
small, certainly under 5 percent, preferably under 1 percent. It is good prac- 
tice to check the quality of the discretization by observing or computing the 
degree of continuity in the stress field. 

2. Compute the equivalent elastic strain f* in each Gauss point and let (f)* 1 ) = 

3. Using (?)<*), compute the secant modulus corresponding to each Gauss 
point from the one-dimensional stress-strain curve. 

4. In each Gauss point for which i > t Y determine the elastic-plastic material 
stiffness matrix. Recompute the stiffness matrices for those elements for which 

i > «y in one or more Gauss points, and obtain a new finite element solution 

Ak+l) 

U FB • 

5. Using 25i fc * and l *, compute the stress tensor components {<r( fc+1 >} in each 

Gauss point, using the total strain computed from and the elastic- 

plastic material stiffness matrix. Determine the elastic strains from {tr( fc+1 >} 
and the elastic part of the material stiffness matrix, i.e., Hooke’s law. Com- 
pute the plastic strain by subtracting the elastic strain components from the 
corresponding total strain components. 

6. Compute equivalent strains (r)( fe+1 ), (?»)( fc+1 ) and s< fc+1 ) from (2.a,b,c). If in 
each Gauss point the following criterion is met: 


- r ' 


( 12 ) 


where r e ia a pre-specified tolerance, then stop, else using compute 

£, (fc+1) , increment Jfc -+ Jfc+ 1 and return to step 3. 
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EXAMPLES 

The solutions of representative examples are solved in the following. Results 
obtained by application of the deformation theory and numerical solution by the 
p-version are compared with results obtained by applications of the incremental 
theory of plasticity and solutions obtained h- and p-extensions. 

The boundary conditions are described in terms of the normal (resp. tangen- 
tial) displacement vector component u n (resp. u«) and the normal (resp. tangential) 
traction vector component T n (resp. T t ). 

Example 1: Plane Strain. 

In this example differences in computed data attributable to alternative 
elastic-plastic models are examined under tight control of the discretization er- 
rors. Two models of elastic-plastic material response are compared: Model 1 is 
based on the deformation theory of plasticity and the von Mises yield criterion, im- 
plemented as described in this paper. The numerical solution was obtained by the 
finite element analysis program PEGASYSf. Model 2 is based on the incremental 
theory of plasticity and the von Mises yield criterion. The numerical solution for 
Model 2 was obtained by an experimental computer program, called FEASIBLE 
[8]. Both programs have p-extension capabilities. 

The solution domain and finite element mesh are shown in Figure 2. Along 
AB and DE symmetry boundary conditions are applied, that is, u„ = Ti = 0. Along 
BC T n = 24, T t = 0. Along CD T n = 30, T t = 0 and along EA T n = T t = 0. 

The material is assumed to be elastic-perfectly plastic. Therefore three pa- 
rameters characterize the stress-strain relationship: The modulus of elasticity (E) 
is 1000, Poisson’s ratio (i/) is 0.3 the yield stress (cry) is 20. The thickness is unity. 

The numerical solutions were obtained by the p-version of the finite element 
method using the six-element mesh shown in Fig. 2 and the product space. The 
product space of degree p is the span of the set of monomials Cn’t * ,j = 0,1,2 ,..., p 
on the standard quadrilateral element (|£| < 1, |r;| < l). For both models the number 
of Gaussian quadrature points was fixed at 14 x 14 for all p-levels. 

f PEGASYS is a trademark of Engineering Software Research and Develop- 
ment, Inc., 7750 Clayton Road, St. Louis, MO 63117. 
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Fig. 2. Solution domain and 6-element mesh for Model 1. 

The number of degrees of freedom (N), the potential energy computed from 
the finite element solutions (il^s) and the estimated relative error in energy norm 
are given in Table 1 for the linear solution. It is seen that the numerical error at 
p = 8 is less than 0.01 percent. 

Table 1. Estimated relative error in energy norm 
for the linear solution. 


p 

N 

Ilfs 

Mm 

_(%)_ 

1 

16 

- 39.9521346350 

7.13 

2 

56 

- 40.1432844703 

1.80 

3 

120 

- 40.1554486036 

0.44 

4 

208 

- 40.1561743285 

0.12 

5 

320 

- 40.1562305887 

0.03 

6 

456 

- 40.1562347043 

0.01 

7 

616 

- 40.1562350709 

0.00 

8 

800 

- 40.1562351098 

0.00 

oo 

oo 

- 40.1562351165 

0.00 


The tolerances for the errors in the iterative solutions were set so small that 
the approximation errors can be considered negligible. Therefore the results show 
variations due to the alternative mathematical models of elastic-plastic material 
behavior. In Model 1 r, = 0.001 was used (see eq. (12)). In Model 2 the error 


- 12 - 


tolerance was set on the residuals in the equilibrium iteration. Specifically, the 


tolerance was y/(Sr) T 6r/r T r < l.o -4 where Sr A = \SK\x — r, [5/C] is the change in 
stiffness matrix after the current iteration, x is the current solution vector and r 
is the current load vector. 


Table 2. Circumferential strain (e t ) 
at the perimeter of the circular hole. Product space, p=8. 


8 

degrees 

Model 1 

Model 2 

Relative 
diff. (%) 

0 

0.1340 

0.1318 

+1.7 

22.5 

0.1272 

0.1248 

+1.9 

45 

0.1029 

0.0989 

+4.0 

07.5 

0.0604 

0.0570 

+4.9 

90 

0.0354 

0.0362 

-2.2 


Models 1 and 2 are compared on the basis of the circumferential strain (e t ) 
along the perimeter of the circular hole. By definition: 


«« 


def 

2 



cos 26 — 



sin 2 8 


where 6 is the angle measured from the positive x-axis. The results are listed 
in Table 2. The results for Model 2 were used for reference, in computing the 

relative differences shown in Table 2. It is seen that the differences between the 
deformation and incremental theories of plasticity are not greater than the errors 
in physical experiments on the basis of which alternative yield criteria are tested 
and the requisite material properties are determined. 


Table 3. p-Convergence of the circumferential strain (« t ) 
at the perimeter of the circular hole. Model 1, product space. 


p N 9 = 0 


2 56 0.1383 

3 120 0.1346 

4 208 0.1329 

5 320 0.1340 


6 = 22.5° 6 = 45° 


0.1112 0.0927 

0.1258 0.1013 

0.1268 0.1004 

0.1269 0.1012 


8 = 67.5° 8 = 90° 


0.0610 0.0407 

0.0593 0.0372 

0.0606 0.0372 

0.0604 0.0366 
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Control of the discretization errors by p-extension was found to be very effec- 
tive. Letting r t = 0.01, and using the default values for the number of quadrature 
points, which are p-dependent, it was possible to obtain substantially the same 
results as in the case of the extremely tight control of discretization and iteration 
errors described above. The results are shown in Table 3. 

A similar problem was solved by Galin using classical methods [7], [8]. Galin 
considered a circular opening in an infinite elastic-perfectly plastic medium, sub- 
jected to uniform stresses at infinity. Galin ’s solution is based on the Tresca yield 
criterion. Because there are differences in the boundary conditions, as compared 
with Models 1 and 2, strict comparison between Galin’s solutions and the solutions 
presented herein is not possible. Nevertheless, because the domain is much larger 
than the circular opening, the differences caused by the differences in boundary 
conditions are very likely to be minor and therefore differences in the solution are 
caused primarily by the differences in the von Mises and Tresca yield conditions. 
The contour lines which separate the yielded and unyielded materials for the Galin 
solution and Models 1 and 2 are shown in Fig. 3. For Galin’s problem this contour 
is an ellipse with major axis of 3.05/r o , 1.64 /r 0 where r 0 is the radius of the circular 
hole. 



Fig. 3. Contours separating the yielded and unyielded regions. 



Example 2: Plane Stress. 

In this problem we consider the elastic-plastic response of a thin perforated 
strip of a strain-hardening material to loading by enforced displacements. The 
results of the numerical analysis are compared with those obtained experimentally 
by Theocaris and Marketos in [11]. 



Fig. 4. Perforated strip. Notation. 

The strip is shown in Fig. 4. All dimensions are in inch units. Taking ad- 
vantage of symmetry, the solution domain was one fourth of the strip which was 
discretized using three finite elements, as shown in Fig. 4. Along AB and DE 
symmetry boundary conditions were prescribed («„ = T t = o); along BC normal 
displacement (a) was imposed (u„ = A, T t = o); along CD and EA the boundary 
was stress free (T n = T« = 0). 

The material properties are typical of an aluminum alloy with yield strength 
in tension c Y = 34, 500 psi, and ultimate strength cuts = 40, 000 psi. The modulus 
of elasticity is E = 9.956 x 10® psi, Poisson’s ratio v = 0.30, and the plastic tangent 
modulus Et = 3.2 x 10 s psi. Plane stress conditions were assumed. 

The stress-strain curve in uniaxial tension, shown in Fig. 5, was characterized 
by five parameters: The slope of the linear part (£?), the slope of the constant 
strain-hardening part (jE<) , the stress at the end of the linear slope (tri = 26,000 
psi), and the smallest values of stress {a 2 = 35,500 psi) and strain (e 2 = 0.0055) 
corresponding to E t . 

For each value of the imposed displacement (a), the resultant force F along the 
edge of the strip and the maximum strain e, at the point of first yield (which occurs 
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Stress (ksi) 



Fig. 5. Example 2: Stress-strain curve in tension based on five parameters, 
at the edge of the hole) were computed from the finite element solution obtained 
for polynomial degree 8 using the trunk space (also known as the ‘serendipity’ 
space). The trunk space of degree p is defined on the standard quadrilateral 
element (|£| < 1, |r/| < l) as the span of the set of monomials *,/ = 0,1,2,..., p, 

i + j < p, augmented by the monomials ^ij, £fj p for p > 2 and by the monomial 
for p = 1. 

The number of degrees of freedom was 211. The estimated relative error 
in energy norm of the starting (linear) solution was 0.23 percent The stopping 
criterion for the nonlinear solution was set at r c = 0.001 (see eq. (12)). 

Table 4. Results for the perforated strip shown in Fig. 4. 


A (in) 

0.0050 

0.0100 

0.0125 

0.0150 

0.0175 

0.0200 

0.0225 

0.0250 

0.0275 

0.0300 


& AV / &Y 

0.217 

0.433 

0.541 

0.645 

0.744 

0.836 

0.917 

0.981 

1.027 

1.055 


Et x f<TY 

0.469 

0.948 

1.251 

1.686 

2.188 

2.659 

3.273 

4.064 

5.269 

6.862 


The results of the analysis for various values of the imposed displacement are 
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presented in Table 4 and in Fig. 6. Figure 6 also includes the experimental results 
which were extracted by reading the values from the plots provided in [10]. The 
normalized stress, a A v fay , is defined as the ratio between the average stress and 
the yield strength: 

&av _ F 
<jy -Amin or- 

and the normalized strain is defined as the ratio between the strain and the 
yield strain (a Y /E). The plastic region is confined up to the maximum normalized 
strain reaching a value of approximately 4.0. 

Normalized Stress 



Fig. 6. Perforated strip: Average stress vs. Maximal normalized strain 

It is seen in Fig. 6 that the computed strain is larger than the experimentally 
observed strain. Other investigators reported similar discrepancies. A possible 
explanation is that in the case of the numerical solution the strain is reported with 
an infinitesimal gauge length wheres experimentally determined strains invariably 
involve some gauge length of finite size. 

Example 3: An axiaymmetrie problem. 

In this problem we consider the elastic-plastic behaviour of a thin-walled 
spherical pressure vessel with a cylindrical nozzle under uniform internal pressure. 
The results of the analysis for an elastic - perfectly plastic material are compared 
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with those obtained experimentally by Dinno and Gill in [12], and numerically by 
Zienkiewicz in [13]. 

The generating section and the finite element mesh, consisting of 14 elements, 
are shown in Fig. 7. All dimensions are in inch units. The material properties 
are typical of a steel alloy with yield strength in tension oy = 40, 540 psi, modulus 
of elasticity E = 29.12 x 10® psi, Poisson’s ratio v = 0.30, and zero strain hardening 
{Et= 0 ). 



Uniform pressure (T n = -p, T t = o) was imposed on the inner surface of the 
vessel. The external surface was stress free (T„ = T t = o). The displacement 
constraints (u r = u, = o) are indicated in Fig. 7. 

The objectives of the analysis were to determine the vertical displacement 
of point A (u?) for a range of pressure values which cause the vessel to yield 
extensively and to determine the size and shape of the resulting plastic zone. 

A sequence of linear solutions was obtained by p-extension using the trunk 
space. The estimated relative error in energy norm of the finite element solution 
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Table 5. Results for the pressure vessel (Example 3). 


Pressure 

Diapl. u? 

(pai) 

(in) 

760 

7.02 x 10~ 3 

900 

8.51 x 10" 3 

1000 

10.78 X 10~ 3 

1080 

15.34 x 10 -3 

1120 

19.66 x lO" 3 

1140 

22.45 X 10 -3 

1160 

27.91 X 10 -3 

1180 

35.65 x 10 -3 

1200 

51.44 X 10 -3 


1400 
1200 
1000 
800 
600 
400 
200 
0 

0.0 10.0 20.0 30.0 40.0 50.0 60.0 

Vertical Deflection of A (x 1000 in) 

Fig. 8. Example 3: Internal pressure vs. axial displacement uf. 

at polynomial degree of 8 (trunk space), was 1.0 percent. There were 1056 degrees 
of freedom. The nonlinear analysis was performed at p-level 8 with the stopping 
criterion r e = 0.01 (see eq. (12)). The results of the analysis for various values of 
the internal pressure are presented in Table 5 and in Fig. 8 which also includes the 
experimental and the finite elements results given in [12]. It is worth noting that 
the plastic zone spreads over the entire section of the nozzle-sphere intersection for 
values of p > 900 psi (see Fig. 9). Good agreement with the experimental results 
was obtained even for high values of pressure. The boundaries of the plastic zone 
for various values of the applied pressure are shown in Fig. 9. These results are 
substantially the same as those presented in reference [13]. 



- 19 - 





Fig. 9. Example 3: Boundaries of the plastic zone for various pressure values. 
Example 4: Limit load in the case of plane strain. 

In reference [5] Nagtegaal, Parks and Rice observed that finite element so- 
lutions based on the displacement formulation exhibit much too stiff response in 
the fully plastic range. Consequently, finite element solutions often exceed the 
limit load by substantial amounts and in some cases have no l imi t load at all. 
This is because plastic deformation occurs at a constant volume. In h -extensions 
based on the displacement formulation the constant volume constraints grow at 
the same or comparable rate as the number of degrees of freedom, hence locking 
occurs. This point is discussed in some detail in [6] also. Locking does not occur in 
p-extensions, however [6], [14]. The following example demonstrates that the for- 
mulation described in this paper will give the correct limit load when p-extension 
is used. 

The most challenging example presented in reference [5] is the computation of 
the limit load for a deep double-edge-notch (DEN) plane strain tensile specimen. 
The example is challenging because the crack is very deep, the ligament is only 
1/9 th of the crack size, hence the crack tip singularity is strong. 

The solution domain is shown in Fig. 10. The boundary conditions are as 
follows: On segments AB and EA symmetry conditions are prescribed (u„ = T t = 
o); on segment BC uniform normal displacement is imposed (u n = S/2, T t = o); 
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segments CD and DE are stress free (T„ = T t = o). The modulus of elasticity and 
Poisson’s ratio were 1.0 and 0.3 respectively, the yield stress was 1.0 also. 

The finite element mesh was graded in geometric progression toward the sin- 
gular point using the grading factor of 0.15. In addition, the mesh quality was 
checked by plotting the stress contours in the course of the nonlinear computa- 
tions. The smoothness of the stress contours across interelement boundaries is an 
indicator of mesh quality. With the exception of the final load case, the trunk 
space at p = 8 was used. The number of degrees of freedom was 2143. In the 
final load case the product space was used at p = 8. The corresponding number of 
degrees of freedom was 4183. 


Net Stress/Yield Strength 



Fig. 11. Example 4. Normalized average stress vs. normalized displacement. 

For the von Mises yield condition the limit load in terms of the net stress on 
the ligament is given by <r« m = (2 + jr)oy/\/3 « 2.97oy . The normalized stress (i.e., 
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the average stress on the ligament divided by ay is plotted against the normalized 
displacement in Fig. 11. The normalized displacement is defined by E6/(ayw) 
where w is the width of the strip (in this example tv = 2.0). It is seen that the 
theoretical limit load is reached at high values of the imposed displacement. 


- 22 - 



SUMMARY AND CONCLUSIONS 

The deformation theory of plasticity is a useful model of elastic-plastic be- 
havior under certain restrictive assumptions which were described in this paper. 

On comparing the deformation theory with the incremental theory in a care- 
fully controlled numerical experiment, Example 1, it was found that the differences 
are smaller than errors in physical experiments designed to test alternative hy- 
potheses concerning elastic-plastic constitutive relationships and the natural vari- 
ations in elastic-plastic material properties. In this experiment the plastic zone 
was completely confined by an elastic zone and the exact solution was smooth. 

The other examples point to the same conclusion: In the numerical exper- 
iments described under Examples 2 to 4 the accuracy of the reference solutions 
was not known, nevertheless it can be said that the results obtained by means 
of the deformation theory were within a few percent of results reported by other 
investigators based on the incremental theory. 

While of course general conclusions from a few experiments would not be 
warranted, the results are consistent with those obtained by Hodge and White [l] 
and with Budiansky’s observation that “deformation theories of plasticity may be 
used for a range of loading paths other than proportional loading without violation 
of the general soundness of a plasticity theory” [2]. In fact, the differences are so 
small that physical experiments could not distinguish between the deformation 
theory and the incremental theories of plasticity in any of the examples discussed 
in this paper. 

The formulation of the elastic-plastic problem described in this paper is based 
on the displacement method. This is possible because p-extensions are not sus- 
ceptible to Poisson ratio locking. The correct limi t loads are obtained. 

It has been demonstrated for the deformation theory of plasticity that p- 
extensions are effective for controlling errors of discretization associated with 
elastic-plastic material behavior. Similar results have been obtained for the in- 
cremental theory even in cases where unloading resulted in reverse plasticity [15], 
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